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ABSTRACT’ 


Higher-order modal methods for predicting thermal and 
structural response are evaluated. More accurate methods or ones 
which can significantly reduce the size of complex, transient thermal 
and structural problems are desirable for analysis and are required 
for synthesis of real structures subjected to thermal and mechanical 
loading. A unified method is presented for deriving successively 
higher-order modal solutions related to previously-developed, 
lower-order methods such as the mode-displacement and mode- 
acceleration methods. A new method, called the force-derivative 
method, is used to obtain higher-order modal solutions for both 
uncoupled (proportionally-damped) structural problems as well as 
thermal problems and coupled (non-proportionally damped) 
structural problems. The new method is called the force-derivative 
method because, analogous to the mode-acceleration method, it 
produces a term that depends on the forcing function and additional 
terms that depend on the time derivatives of the forcing function. 


The accuracy and convergence history of various modal methods 
are compared for several example problems, both structural and 
thermal. The example problems include the case of proportional 
damping for: a cantilevered beam subjected to a quintic time- 

varying tip load and a unit step tip load and a multispan beam 
subjected to both uniform and discrete quintic time-varying loads. 

Examples of non-proportional damping include a simple two-degree- 
of- freedom spring-mass system with discrete viscous dampers 
subjected to a sinusoidally varying load and a multispan beam with 
discrete viscous dampers subjected to a uniform, quintic time- 
varying load. The last example studied is a transient thermal 
problem of a rod subjected to a linearly-varying, tip heat load. 

The higher-order modal methods are shown to converge to an 

accurate response using fewer eigenmodes than lower-order modal 
methods. The force-derivative method is very effective in 

representing the response of the important, but otherwise neglected, 
higher modes for structural problems in which there are a large 
number of closely-spaced frequencies (e.g., a multispan beam or a 
large area truss-type structure). However, for response times close 
to discontinuities in the forcing function and/or its derivatives, the 
mode-acceleration or force-derivative methods must include 
appropriate jump conditions or additional modes to insure accuracy. 
The higher-order modal methods are also effective in solving 
transient thermal problems efficiently. 
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Chapter 1 


Introduction 


Transient thermal and structural analysis of complicated 
engineering problems which are modeled as discrete multi-degree-of 
freedom systems often require the solution of very large systems of 

equations. The accurate solution of very large systems of coupled 

differential equations may be difficult and computationally 
expensive. In addition to calculating the transient response, 
optimization of the structure under transient thermal and mechanical 
loading may be necessary. Automated structural design requires 
numerous solutions of the coupled system of equations to assess the 
effects of changes in design parameters on response quantities and 
constraint boundaries. Since the problem is transient, continuous 
constraint boundaries in time have to be approximated, often 

requiring many discrete times to be used to approximate the 

constraint boundaries and ensure that the critical times when the 
maximum response occurs are properly represented. Automated 


design of even moderate-size problems can become intractable when 
designing for transient loading conditions. Hence, reducing the size 
of such systems is highly desirable from the standpoint of increased 
computational efficiency. Some of the many methods for reducing 
the size of discrete multi-degree-of-freedom structural dynamic 
systems include mass condensation methods (e.g., refs. 1 to 7) and 
reduced basis methods (e.g., refs. 8 to 15). 

The mass condensation methods typically reduce the size of the 
dynamic system of equations by condensing out less-important 
degrees of freedom from the problem. One of the most popular 
methods, proposed by Guyan (ref. 1), assumes the relationship 
between the dependent and independent degrees of freedom can be 
determined by the static relationship between them. Other mass 
condensation methods either approximate a relationship between 
dependent and independent degrees of freedom of the dynamic 
system (ref. 7) or solve exactly for this dynamic relationship (ref. 3). 
To achieve reasonably accurate results, the set of dependent 
degrees-of-freedom must be chosen carefully or some of the lower 
frequencies in the eigenspectrum may be lost. It was also shown in 
reference 7 that large errors can occur when Guyan reduction is used 
to expand the mode shapes and calculate forces. Modified mass- 
condensation methods have been developed to help reduce these 
deficiencies (refs. 3, 6, and 7). An attempt to extend the Guyan 
reduction method to transient thermal analysis (ref. 16) resulted in 
the conclusion that a crude thermal finite-element model (having 
fewer degrees of freedom than the model reduced using Guyan 


reduction) produced results as good as or, in most cases, better than 
those obtained using Guyan reduction. 

The reduced basis methods, on the other hand, use an expansion 
theorem (ref. 17, pp. 90-91) and a reduced set or subset of linearly 
independent basis vectors to approximate the transient response. 
These methods use either a truncated set of complete basis vectors 
(e.g., eigenmodes (ref. 18), Ritz vectors (ref. 8), or Lanczos vectors 
(ref. 19)) or a combination of basis vectors (e.g., eigenmodes and Ritz 
vectors (ref. 14). When a reduced basis method is based on the 
eigenmodes of the system, the method is referred to as a modal 
method. One of the first reported uses of eigenmodes to reduce the 
size of a problem and calculate the dynamic response of an elastic 
structure was by Biot and Bisplinghoff (ref. 18). In that classic paper, 
the authors recognized that the transient, linear response of 
structures can be calculated by a superposition of the natural 
oscillations or modes since the eigenmodes are orthogonal and 
linearly independent and form a complete set (this method is herein 
referred to as the mode-displacement method (MDM)). In addition, 
for linear dynamic problems without damping or for proportionally- 
damped problems, the equations of motion of a multi-degree-of- 
freedom system uncouple and can be solved individually as single- 
degree-of-freedom systems. 

The usefulness of any reduced-basis or mode-superposition 
method lies in its ability to predict accurately the transient response 
using only a small number of basis vectors or modes. The need for 
only a small number of basis vectors for an accurate response is 



especially well-suited for modal methods, where it is 
computationally expensive to calculate the eigenvectors or modes of 
large systems as compared to other, less expensive reduced basis 
methods, which use Lanczos vectors (ref. 19) or Ritz vectors (ref. 
14). Toward this end, the MDM has proven useful in calculating 

transient structural displacements for most dynamic structural 
problems, using only a small percentage of the lower frequency 
mode shapes. 

The MDM can be considered to be a generalization of a Fourier 
series approximation. It is well known that Fourier series 
representation of discontinuous functions converge slowly, and that 
the derivatives of the series may not converge at all. This type of 
convergence problem is known as the Gibbs phenomenon (ref. 20). 
Hence, the MDM can be expected to converge slowly when the 
applied loads exhibit discontinuities in time or space. This slow 
convergence is exhibited in the transient response problem of a 
string with a point load (ref. 21), and other examples (e.g., refs. 22 
and 23). The MDM is also less accurate in predicting transient 

stresses (refs. 9 and 24) and, hence, requires more modes to obtain 
converged accurate stresses. This decreased accuracy in predicting 
stresses is understandable since the stresses are related to the 
spatial derivatives of the displacement vector and, hence, errors in 
the displacements become magnified upon differentiation. In 
addition, some structural problems which have closely-spaced 
natural frequencies, such as large area space structures or multispan 
beams, experience very slow convergence and require a large 



number of modes to predict even the displacement response 
accurately. 

The mode-acceleration method (MAM) was developed as a means 
to improve the convergence of the MDM. The origin of the method is 
attributed to Williams (ref. 25), but as pointed out by de Veubeke 
(ref. 9), the original concept was stated by Lord Rayleigh (ref. 26) as 
early as 1877. Lord Rayleigh noticed that when the periods of the 
forces operating on a system are long relative to the free vibration 
periods of the system, the inertial forces of the system can be 
neglected. The MAM was popularized by Craig (ref. 27) in 1981 and 
put into a more familiar form. The MAM improves the low- 
frequency or pseudo-static response convergence because it 
incorporates, as a separate term, the pseudo-static response in the 
solution which approximates, to some degree, the flexibility of the 
higher modes which are neglected in the modal summation. Maddox 
(ref. 24) indicates an improved accuracy in dynamic force 
calculations using a MAM-type approach as compared to the MDM 
approach. Hansteen and Bell (ref. 28) indicate that the inaccuracies 
of modal truncation in the MDM can be caused by components of the 
load which are orthogonal to the eigenmodes included in the solution 
and show that Maddox's proposal is equivalent to the inclusion of an 
approximate expression for these components of the load. 
Anagnostopoulos (ref. 29) indicates deficiencies in both the MDM and 
MAM in calculating forces and stresses in offshore structures when 
subjected to earthquake-type loading. Leger and Wilson (ref. 30) 
compare the numerical efficiency of several forms of the MAM as 


presented in references 28 and 29. An in-depth comparison of the 
MDM and MAM methods (ref. 11) indicates that the MAM converges 
to an accurate solution with fewer modes than the MDM. This study 
includes a comparison of the effect of proportional damping level and 
load frequency on the relative accuracy of the response using the 
MDM and MAM methods. 

Attempts at using mode-superposition methods to solve transient, 
linear thermal problems (refs. 31 to 35) have been unsuccessful 
because thermal problems exhibit a wide spectrum response where 
very high frequencies are excited and, hence, a prohibitively large 
number of "thermal modes" are necessary for an accurate solution. 
The use of "thermal modes" was first suggested as a solution to 
transient thermal problems by Biot in reference 36. In that paper, 
Biot develops an analogous MDM solution and also includes what is 
analogous to a MAM solution for a plate with one edge insulated and 
the other edge suddenly brought to constant temperature. 

The work of Ramberg (ref. 37) reveals that the MAM can be 
derived from the MDM by integrating the convolution integral 
portion of the solution by parts once with respect to time. Leung 
(ref. 12) improves the convergence of the MAM for undamped 
systems by integrating the convolution integral several more times. 
A similar method for developing highly convergent modal solutions 
to linear dynamic structural problems was suggested by Likhoded 
(ref. 38), which the author calls "multiple extractions of the pseudo- 
static component". Likhoded presents a recursive technique for 
developing higher-order or faster-convergent modal solutions. 


Borino and Muscolino (ref. 39) developed a modal method, called the 
dynamic-correction method (DCM), which separates the 
homogenuous and complementary solutions of the transient 
response. The DCM assumes that the transient response at time, t, 
can be considered as the sum of a pseudo-static response or 

particular solution, which depends on the load at time, t, plus a 
dynamic correction which takes into account the dynamic solution 
due to the loads from the initial time to time t. 

The present study extends the work of Leung (ref. 12) by 
developing a unified method (refs. 13 and 40) for deriving higher- 

order modal methods which are general in that they: (1) are easy to 

implement into existing computer programs, (2) can represent 
proportional as well as non-proportional damping, and (3) can be 
extended to solve transient, linear heat transfer problems. The 
newly developed method is called the force-derivative method (FDM) 
because, analogous to the MAM, the FDM produces a term which 
depends on the forcing function and additional terms which depend 
on the time-derivatives of the forcing function (ref. 13). These 
additional terms produce successively higher-order approximations 
to the contributions of the higher modes which are neglected in the 
modal summation. 

The effects of various factors on the rate of convergence or 

accuracy of the various modal summation methods is investigated. 
These factors include: the differentiability of the forcing function, the 
frequency of the forcing function, the level of damping, and the time 
at which the response is calculated. The forcing functions were 


chosen to illustrate the effect of a temporally continuous forcing 
function which has continuous higher derivatives, a discontinuous 
forcing function (a unit step function), a spatially continuous loading 
distribution, and a spatially discontinuous loading. Structural 
problems include both proportional and non -proportional damping as 
well as the effect of discrete dampers. A quintic function of time was 
selected to illustrate what happens when the force or one of its 
derivatives vanish at some point in time. The convergence for a 
sinusoidal forcing function is also studied. 

A series of numerical examples has been selected to illustrate the 
adequacy and/or inadequacy of each method. The examples, which 
assume proportional damping, are: (1) a uniform cantilevered beam 

subject to a quintic time-varying tip load, and a unit step tip load 
and (2) a multispan beam (10 equal-length spans) subjected to both 
uniform and discretely applied, quintic time varying loads. Examples 
assuming non-proportional damping include: (1) a simple two- 

degree-of-freedom spring-mass system with discrete viscous 
dampers subjected to a sinusoidally varying load on one of the 
masses, and (2) a multispan beam with discrete viscous dampers 
subjected to a discrete, linearly-varying load. In addition, the 
suitability of using higher-order modal methods to solve transient, 
linear heat transfer problems is investigated. The last problem 
studied is a rod heated at one end by either a step tip heat load or a 
linear time-varying heat load while the other end is maintained at a 
uniform temperature. 



A method for developing higher-order modal methods which 
mathematically unifies previously developed improved modal 
methods (refs. 10 to 12, 38, and 39) is presented in Chapter 2. 
Methods for solving proportionally and non-proportionally damped 
systems of equations using both damped and undamped eigenmodes 
are presented. The first-order, damped-mode formulation of the 
equations of motion is especially suited for solving transient heat 
transfer problems. Two different error norms are used to evaluate 
the various methods and are presented in Chapter 3. Analytical 
results are divided into two separate chapters, one for Structural 
Analysis (Chapter 4) and another for Thermal Analysis (Chapter 5). 
A summary and conclusions are given in Chapter 6. 


Chapter 2 


Unified Derivation of Modal Methods 


As mentioned in the introduction, there have been many modal 
summation methods and forms of implementation of such methods 
presented in the literature (e.g., refs. 3, 9, 10-12, 14, 18, 24, 25, 28, 
29-31, 36, 38, and 39). The present chapter will describe a unified 
method for deriving the simplest or zeroth-order form of the modal 
summation methods, the mode-displacement method (MDM). The 
method will be shown to be useful for deriving higher-order 
methods such as the mode-acceleration method (MAM) (which will 
be shown to be a first-order method). This unified method, which is 
herein called the force-derivative method (FDM), can be developed in 
two forms, a first-order or damped-mode formulation which 
operates on the equations of motion in first-order form and a 
second-order or natural-mode formulation which operates on the 
equations of motion which are in a second-order form. The former 


uses the damped system eigenvalues and eigenvectors, whereas the 
latter uses the natural modes and eigenvalues of the system. In 
addition, alternate forms of the above-mentioned formulations are 
presented which are easier to implement into existing computer 
programs. 


2.1 First-Order or Damped-Mode Formulation 


The equations of motion, in matrix form, of an n-degree-of- 
freedom system, together with the initial conditions, are given by 


Mii + Cu + Ku = Q(t) (1) 

u(0) = u 0> u(0) = u 0 


where M, C, and K are the mass, damping, and stiffness matrices of 
the system; u and Q are the displacement and load vectors, 
respectively, and a dot denotes differentiation with respect to time. 


Transforming eq. (1) into first-order form results in the following 
system of equations: 


MY + K Y = Q Y(0) = Y o (2) 


where 



Assuming a solution to the homogeneous form of eq. (2) as: 

a r t 

Y(t) = e <D r ( 3 ) 


results in an eigenvalue problem 

a r M <D r + K 4> r = 0. (4) 


12 


The eigenvectors (<I> r ) are normalized such that 
T _ 

d> r M <D r = 1.0 

and then (5) 

T _ 

O r K <J> r = -a r 

Equation (5) can also be written in matrix form as 


[<D] T M [O] = fl] and [0] T K [O] = [a] 


where [®] is a 2n-by-2n modal matrix with its ith column equal to 
O r and [a] is a diagonal matrix consisting of the a r 's. 


A solution to eq. (2) is assumed in the form of the following modal 
summation 

2 ii 

Y(t) = X^rZrCt) (6) 

r= 1 


1 3 


T 

Substituting eq. (6) into eq. (2) with premultiplication by O r 
results in the following, uncoupled, system of equations: 


Z r ■ Z r - Q 


Z r (0) = Z ro = 4> r M Y 0 


( 7 ) 


The solution to eq. (7) is 


a r t f a r (t-x) T _ 

Z r (t) = Z ro e + J e 0> r Q(x)dx ( 8 ) 

o 


Hence, the solution of eq. (2) becomes 


Y(t) = 



a r t 

e 


r= 1 


l 

f Otr(t-X) 

+ Je F r (x)dx 


where 


Fr(x) = d> T r Q(x) 


( 9 ) 


If the forcing function has continuous derivatives, the convolution 
integral of eq. (9) can be integrated by parts to produce higher-order 

1 4 


modal methods (ref. 13). For example, if it is integrated by parts 
once, the following expression results 



( 10 ) 


If all the modes are used in the second-to-last term in eq. (10), 
this term can be written as 

2 n ^ 

X ^r—PrO) = ['!>] [<J>] T 0(t) = K Q(t) (11a) 

r= 1 

since from eq. (5) 

[[«*>! K [<&]] =-[«] or [4>] R [O ] T = [^] (lib) 

T 

Pre- and post-multiplying eq. (lib) by [*1 and [* . 
respectively, leads to 

K = -fO] ^ [<1)] T (11c) 

1 1 1 

where ~ is a diagonal matrix whose elements are . 

(XJ tXj- 

1 5 



If equation (11) is substituted into equation (10), the resulting 
expression can be written as 


2 n 


s 


Y(0= / , O r 


F r (oy 


Zro + a 


r J 


a 


t O r f a r (t-x) 
+ 77 " Je 


F r (x)dx 


r= 1 


_-l_ 

+ K Q(t) 


( 12 ) 


If the forcing function has continuous derivatives up to order N-l 
(C N '' continuity), the convolution integral of eq. (9) can be integrated 
by parts N times, resulting in the following expression: 



i=l 


where the superscript (i-1) denotes the (i-l)th derivative with 
respect to time. 


If all the modes are used in the last N-terms of eq. (13), these 
terms can be represented as functions of M and K (ref. 40) as follows 



i=l 


2.2 Alternate Damped- Mode Formulation 

Equations 8 to 14 are analogous to those presented in reference 
13. However, the present expressions solve a first-order system of 
equations, using the damped modes, 0 r , to decouple a non- 
proportionally damped system (the damping matrix is not a linear 
function of the mass and stiffness matrices). The first-order system 
of equations is twice as large as the second-order system. Equation 
14 represents a means for developing higher-order modal methods 
than either the MDM or MAM. This method is called the force- 
derivative method (FDM) (refs. 13 and 40) because it produces terms 
which are related to the forcing function and its time derivatives. 



Equation 2 can also be considered to represent a heat conduction 
problem where, for that problem, M represents the capacitance 

matrix, K represents the conductance matrix, Q represents the thermal 
load vector, and Y is the vector of nodal temperatures. 

The MDM uses a subset, m (m < 2n), of the eigenmodes to reduce 
the size of the problem and solves for Z r using eq. (8) (or simply by 

numerically integrating eq. (7)) and substitutes these values into a 

reduced modal summation in eq. (6) to approximate the response, 
Y(t). The MDM can be classified as a zeroth-order method because it 
is equivalent to using the FDM (eq. (14)) with N = 0. An analogous 
form of the MAM uses eq. (12), and a reduced modal summation to 
approximate Y(t) and can be classified as a first-order method 
(N = 1 in eq. (14)). The FDM uses eq. (14) with N > 1. (Results 

presented in this study assume N = 4 when referring to the FDM 

method). Reference 13 showed that the expressions obtained using 

the FDM offer improved approximations to the contributions of 

higher, neglected, modes for several structural problems. 

The FDM (eq. (14)) can be derived using an approach similar to 
that used in reference 39 which results in a form which is more 

suitable for inclusion into existing thermal and structural analysis 

codes. A numerical approach can be derived, similar to that 

presented in reference 39, which approximates the forcing function 

as a piecewise differentiable polynomial and which numerically 
integrates the reduced system of equations (eq. (8)). For example, 



assuming the forcing function is continuous; eq. (7) could be 
rearranged as shown 


Zr(t) = ~ 4>r Q(0 + Z r (t) (16) 

Using eqs. (6), (11), and (16), the response can be written as 

2 n 

Y(t)=K Q(0 + ^® r -J^Z r (t) (17) 

r=l 

The last term can be evaluated using eq. (8) and Leibnitz's rule for 
differentiation of an integral to produce the following: 

t 

T _ a r t f a r (t-x) T _ 

Z r (t)=O r Q(t) + a r Z ro e + a r Je O r Q(x)dx (18) 

o 

Y(t) can be approximated using only a subset of the modes for the 
last term in eq. (17), and using eqs. (8) and (18), eq. (17) becomes 


1 9 


Y(t) s (k + <t> & ’o T )q( 0 + $Z(t) 


(19) 



where the A denotes a reduced set of modes (m < 2n), a represents 

a diagonal matrix consisting of a reduced number of eigenvalues 

a/ 1 , and Z(t) can be calculated by either numerically integrating 

eq. (7) or by using an analytic expression, when applicable, for a 

given forcing function. Equation 19 is an alternate form of the MAM. 

If the forcing function can be assumed to be differentiable, 
eq. (7) can be differentiated once with respect to time and re- 
arranged as shown below: 

z r d) = ■ <t> r Q(t) + z r (t) (20) 

Re-arranging eq. (7) and substituting for Z r from eq. (20) results 
in 


Zr(t) — ' " 9> r Q(t) ' 2 ^r Q(0 + 2 (21) 

r rv rt 


Using eqs. (5) and (6), the response can be written as 


1 _ 

Y(t) = K Q(t) - K MK Q(t) + 



( 22 ) 


where the second term in eq. (21) is determined using eq. (11) and it 
is assumed that the modes can be normalized as shown by eq. (5). 
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The last term in eq. (22) can be evaluated using eq. (8) and 
Leibnitz’s rule for differentiation of an integral to produce the 
following: 


T _ Tj. 

Z r (t) = ot r O r Q(t) + 0fQ(O + 


t 


2 art 2 f ar(t-x) T _ 
a r Z ro e + a r Je 0 r Q(x)dx 
0 


(23) 


Once again, Y(t) can be approximated by using a subset of the 
modes for the last term in eq. (22), and using eqs. (8) and (23), eq. 
(22) becomes 


(-• 1 a a -IaT) a a -2aT ) • 

Y(t) = (^K + O a O jQ(t) + (-K MK +0a<J> J Q(t) 

+ ^Z(t) (24) 


2 1 


If a C“ differentiable forcing function is assumed, eq. (7) could be 
differentiated twice and a procedure, similar to that outlined by 
eqs. (20) to (24), would produce the following expression: 

(- X A a -1AT >| ( A a -2aT \ • 

Y(t) = {K + O a d> jQ(t) + 1 : KMK +<Da <D J Q(t) 

f- 1 — 1 — 1 £ a-3 AT >, •• A A 

+ (K MK MK +Oa jQ(t) + OZ(t) (25) 

This expression for Y(t) shown by eq. (25) can be expanded, by 
assuming the forcing function is C N differentiable, to give 

N+I 

A A-iAT|_(i-l>] AA 

[^-K Mj K + a> a <I> jQ(t) J + <SZ(t) (26) 

i=l 

Compared to eq. (14), the alternate formulation of the FDM (eq. 
(26)) does not require the solution of a convolution-type integral. In 
addition, the last term of eq. (26) is identical to a mode-displacement 
solution and, as such, the form of the FDM as given by eq. (26) is 
more suited for inclusion into existing computer codes. 
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2.3 Second-Order or Natural-Mode Formulation 


Expressions which use the undamped natural modes can be 
developed in an analogous manner as the damped-mode formulation 
and result in an expression similar to eqs. (19) and (24) to (26). 
Beginning with the second-order system of equations (eq. (1)), the 
undamped or natural modes, (j) r , are determined by solving the 

following eigenvalue problem: 

K(J> r = (Of M(|) r (27) 

where (D r is the rth circular natural frequency. 

The modes are normalized as follows: 


T 

(t) r M <t) r = 1.0 


so that 


1 2 
(j) r K (j)|* — (Or 

Hence, the displacement response can be represented as 


(28) 


u(t) = X ♦rflrCO ( 2 ^) 

r= 1 
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Using eq. (29) and premultiplying eq. (1) 


by <|> r results in 


2 T 

q + A q + Q q -[♦] Q(t) 


( 30 ) 


where 

A = [<f] Tc K] 

and 

= [4>] T k|>] 

where | <> | is the matrix of undamped eigenmodes, Q is a diagonal 

2 2 

matrix whose diagonal terms can be represented as £2 ii = G)j and, for 
proportional damping (where the damping matrix can be 
represented as a linear combination of the mass and stiffness 
matrices), A is also a diagonal matrix whose diagonal terms can be 
represented as Ajj = 2£j(0j. 

Assuming proportional damping and zero initial conditions, the 
solution to eq. (30) can be written as 


i 


q 



CrCOr(t-'C) 

e sin 


T 

WdrU-'0<l ) rQ('t) dt 


(31) 


where 


®dr = 



(Cr“>r ) 2 
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The MAM is equivalent to the FDM of order one (assuming one 
integration by parts of the convolution integral). This equivalence 
can be shown by integrating eq. (31) by parts once with respect to 
time, premultiplying by (j) r , and substituting into eq. (29). This 

results in: 


m 

u(t) 5 X <Mr(0 
r=1 


m 


X 


, . -<t>r Q(°) -CrCOrt fCr®r . r . . r ,1 

9r| 2 e 17^7" sin(o)drt) + cos(cod r O| 


CO, 


w d r 


r=l 


+“7 < t> r Q(0) 
(Or 


(32) 
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-jrQjj -CrWrd-^fCrWr 


j 

0 


CO, 


l“d7 Sin<ddr< ‘ t ’ X ^ + COS ®dr(l-'l ; )]cl'C 
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If all the modes are used in the second-to-Iast term of eq. (32), it 
can be written as 


m 


u(t) = K Q(t) + 7 < f) 


I' 


-$ T r Q(0) -Cr© r t fC r O) r . . ^ , . .1 

2 e 103^7 s,n(Q>drt) + cos(codrt)j 


r=l 


(33) 


Q( x ) ■Cr®r( t_ ' c ) fCr ro r . ] 

? 6 s intO dr (t-T) + COSO) dr (t-T)|dT 


(0, 


Equation (33) can be shown to be equivalent to the MAM (ref. 11), 
the expression of which is shown below: 


-1 

u(0 s K Q(t) 



2Cr - 1 •• 

q(t) + — q(t)} 
(0 r 


(0 


(34) 
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Using Leibnitz’s rule for differentiation of an integral and 
differentiating eq. (31) with respect to time twice gives 


1 “Cr®r(UX) « f 

q r (0 = I e -Cr« r sin(0 dr (t-x)+ (O dr cos(O dr (t-x)}(t) r Q(x) dx 


(0 d 


r ^ 

r o 


and 


(35) 


T if ^ r ®r(t-X) w 

q r (0 = <t>rQ(0 + 7T” J e x 

^dr 0 


T 

{(2(C r (0 r ) 2 -(0 2 r )sinw dr (t-x) + 2C r co r (O dr cosa> dr (t-x) }<|> r Q(x)dx 


Substituting eqs. (35) into eq. (34) and simplifying, results in 


m 


-1 

u(t) = K Q(t) 4- 


I 


r=l 


1 t '' 

— ^r 2^r CKO 4* t^r^rO) 

co r 


(36) 


Integrating the convolution integral expression for q r (t) with 
respect to time once, as was done earlier (eq. (32)), it is easy to 
verify that the MAM (eq. (36)) is equivalent to the FDM of order one 
(eq. (33)). This equivalence was also shown in references 13 and 37. 


2.4 Alternate Natural-Mode Formulation 


An alternate natural-mode formulation is developed in a similar 
manner as the alternate formulation of the first-order or damped- 
mode formulation of Section 2.2. This alternate form is well-suited 
for inclusion into existing computer codes. Assuming the forcing 
function is differentiable, eq. (30) can be differentiated twice and 
substituted back into eq. (30) to produce the following expression: 


q(t) = Q 2 [<|)] T Q(t) - Q 2 AQ 2 [(j)] T Q+ [q 2 AQ 2 AQ 2 - Q* 2 q' 2 J<}> T Q 


-2 -2 - 2 -2 -2 -2 -2 1 (3) 

+ Q £2 A - Q AQ AQ A + Q AQ q (t) 


-2 -2 -2 -2 -2l (4) 

Q Q - Q AQ AQ J q (t) 


(37) 


If eq. (31) is differentiated with respect to t and the expressions 
for q(3) and q(4) are substituted into eq. (37), the entire 
expression reduces to the following: 


u(t) 2 (k 
-1 * 1 


-I AA-2aT>| (- 1 -1 A a-2 a a-2aT \ • 

- <1> Q ({) J Q(t) - ^K CK -<J)QAQ(j) Jc 


Q(t) 


K MK 1 - K CK CK ) - ($ & ^ 2 <j> T - $ tfk k Q(t) 


A A 

+ <})q(t) 


(38) 


Equation (38) agrees with results presented in reference 39 and, 
as shown in reference 39, is also valid for non-propordonally 
damped structural systems. Assuming higher-order piecewise 
differentiable forcing functions, the method above can produce 
successively higher-order modal methods and, as such, is just 
another formulation of the FDM. 
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The expression for an N lh -order modal method can be expressed 
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2.5 The Dynamic-Correction Method (DCM) 


The dynamic-correction method (DCM), reference 39, assumes a 
solution to eq. (2) in the form 


Y(t) = Y P (t) + Y c (t) 


(40) 


where Y P (t) is a particular solution of eq. (2) and Y c (t) is the 
complementary solution which represents the effects of initial 
conditions. In modal form, Yp(t) and Yc(t) can be represented as 

Yp(l) = [0] Z P (t) 

and (41) 

Y c (l) = [0] Zc(t) 

where Zp(t) and Z^(t) are the vectors of the particular and 

complementary solutions to the modal coordinate equations (eq. (7)). 

Using eqs. (6) and (41) 

Y(t) = [0] Z(l) = [0][Z(t) - Z P (t)] + Yp(l) 


(42) 



The fundamental principle of the DCM is that, if available, an exact 
particular solution to eqs. (2) and (7) can be used to approximate the 
response (eq. (40)) using a reduced set of modes as shown below: 


Y(t) @ F^tl) + [y P ( t) - F2p(t)] 

(43) 

It can also be shown that in the limit as N goes to infinity, two 
terms in equation (26) can be written as 



Hence, if an infinite number qT integrations-by-parts are assumed 
in the FDM or if the convolution integral vanishes (e.g., for a 
polynomial forcing function of a lower order than the order of the 
FDM) the FDM would be equivalent to the DCM of reference 39. Also, 
if an exact solution to the convolution integrals (eqs. (8), (12), (14), 
(33), etc.) exists and is used, the response can be calculated without 
the errors caused by approximating the forcing function. 
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Chapter 3 


Error Norm Definition 


It is important to develop reliable error estimates to evaluate or 
compare, quantitatively, the various modal reduction methods. As 
mentioned previously, classical modal superposition methods for 
linear systems use a subset of the lower modes to approximate the 
transient response. The error estimates used here compare the exact 
or converged response with an approximation to that response, 
obtained by using a subset of the eigenmodes. These error norms 
are not intended to be used to predict, a priori, the number of modes 
necessary for convergence and, hence, the time for termination of the 
modal series. The convergence of each method (number of modes 
versus the accuracy of the transient response) is measured by using 
one of two relative error norms: a spatial error norm or a time- 

integrated error norm. 


3.1 Spatial Error Norm 


The spatial error norm, e, of an approximation to the temperature 
vector is given by 


e = 



( T - T a ) T ( T - T a ) 
T 

TT 


(45) 


where T represents a converged solution for the temperature vector 
and T a is an approximation which can be based on the first m 
thermal modes. 


3.2 Time-Integrated Error Norm 


A time-integrated error norm, similar to that used in reference 
39, can be used for the error in displacement uj and is shown below: 


J I u j(t) - u?(t)l dl 

£:(%) = X 100 (46) 

1 T 

J I U i( t) I dt 
0 


34 



where £; is the time-integrated error in the ith nodal displacement, 
u j ( t ) is the exact response using all the modes of the system, and 
u?(t) is the approximate response at node i using a subset of the 
lower modes. The time, X, selected as the upper limit of integration 
was chosen to be X = 167t/COf f where w f is the forcing function 
frequency- 



Chapter 4 


Structural Analysis 


The following sections investigate the effect of various forcing functions, 
load distributions, and damping levels on the transient response of a 
uniform cross-section cantilevered beam, a uniform, simply-supported 
multispan beam, and a spring-mass-damper system. The rate of 
convergence of each of the methods (MDM, MAM, FDM, and DCM) is 
expected to depend on the nature of the forcing function, the level of 
damping and the time at which the response is calculated. Proportional 
damping, when assumed, is constant for all modes. The forcing functions 
are selected to investigate the effects of continuous forcing functions with 
vanishing higher derivatives at various times (Q(T) = 1000(T 4 - T^)), the 
effect of a discontinuous forcing function representing a unit step at time T 
= 0 (Q(T) = JLt (T)), a spatially discontinuous forcing function, and the effect 
of a sinusoidal forcing function (Q(T) = sin (OfT), where T is the normalized 

time, T = co () t, and co 0 is the normalizing frequency, to 0 = 

In addition, the effects of non -proportional damping are investigated. It is 
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assumed that for all forcing functions, Q(T) = 0 for T < 0. Several example 
problems, described below, are selected to evaluate the accuracy of each 
method. 

The spatial error norm, e (eq. (45)), is used to evaluate the various 
modal methods. The effectiveness of this error norm in quantifying the 
global error associated with each of the modal methods is demonstrated 
for a cantilevered beam problem. The displacement distribution of a 
cantilevered beam with a quintic varying tip load in time of Q(T) = 
1000(T4 - T^) lb. at a normalized time T = 0.4 and for = 0.05 is shown 

in figure la for each of the various modal methods using only one mode. 
The value of the spatial error, e, is also listed in the figure. Notice that as 
the value of the error norm decreases, the distributions approach that of 
the MDM using 30 modes, as indicated by the solid line. As shown in 
figure la, only one mode is used and for the MDM, e = 0.289 and 
displacement errors are noticeable. However the MAM, FDM and DCM 
results, having error norms of e = 0.0407, 0.0008, and 0.001, respectively, 
are indistinguishable from the converged solution (MDM using 30 modes). 
Similarly, for the normalized moment distribution (Fig. lb) the MDM (using 
the first two modes) has an error e = 0.395, and the MAM, FDM and DCM 
(all three using only the first mode) have errors of e = 0.1190, 0.0023, and 
0.0033, respectively. As shown in figures la and lb, there is a qualitative 
improvement in the solution (response distribution) as the error norm e 
decreases in magnitude. 
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4.1 Proportional Damping 


The equations ot' Chapter 2 are used to study several beam 
example problems: a cantilevered beam with a tip loading and a 
multispan beam with uniform and discrete loadings. These problems 
use analytical expressions for the mode shapes <t> r (x) (refs. 17 and 
41), mode-shape derivatives, and modal coordinates q r (t) (eq. 29). 

The use of analytical solutions to calculate the transient response of 
the modal coordinates, q r , eliminates the need for numerical 

integration in time and associated numerical errors. As mentioned in 
Chapter 2, the MDM and MAM can be considered as the FDM of order 
zero and one, respectively. The results labeled FDM in the figures, 
refer to a fourth -order version of the FDM which is obtained 
assuming the forcing function has derivatives of order four or higher. 
Modal vectors are calculated assuming 51 equally-spaced points 
along the length of the beam. The solution formulation used is the 
second-order or natural-mode formulation (Section 2.3) where the 
convolution integral portion of the expression (e.g., eq. (33) and ref. 
13) is evaluated analytically. 


Moment and shear forces were calculated from the following 
equations: 

M(x,t) = EI-|^ 


and 


(47) 


S(M) = 


For most cases, 30 modes are sufficient for an accurate solution 
and, hence, 30 modes are used to approximate expressions such as 


-l 


K = 



f \ \ 


\®r) 




T 
r > 


-1 -1 
K MK = 



n V T 

<l>r 


4 

v w r y 


, etc. 


(48) 


The multispan beam examples experience much slower convergence 
and, hence, for this problem 50 modes are used to represent the exact 
solution. 


4.1.1 Cantilevered Beam with Tip Loading 


The first problem studied is a uniform cantilevered beam loaded by a 
tip loading. The first 30 natural frequencies of the beam are listed in 
Table 1. The full 30 modes are used in expressions like eq. (48) to 
represent the converged solution. The beam is subjected to various levels 
of modal damping (same value for all modes) and a variety of loading 
conditions. 


4. 1.1.1 Quintic Time-Varying Load.- For this problem the forcing function 
is Q(T) = 1000(T 4 - T 5 ), where T is the normalized value of time. This 
problem was presented in reference 12 for the case of zero proportional 
damping in all modes (Ci = 0). This forcing function was chosen to evaluate 

the various higher-order modal methods since the function or one of its 
derivatives vanishes at various times: Q(T) = 0 when T — 1.0, 

Q(T) = 0 when T = 0.8; Q(T) = 0 when T = 0.6; Q^T) = 0 when T = 0.4; and 

q ( 4 ^T) = 0 when T = 0.2. This fact affects the convergence of the method as 
will be shown subsequently. The forcing function is plotted as a function 
of time, from time T = 0 to time T = 1.2, in figure 2a. The variation of tip 
displacement as a function of time, for the case of zero damping (^j — 0), is 
shown if figure 2b. As shown in figure 2b, the tip displacement decreases 
to a minimum value of about -30.0 in. at 1.1 sec, shortly after the forcing 
function changes sign at T = 1.0 (see fig. 2a). The moment error norm as a 
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function of the number of the modes used in the modal summation is 
shown in figure 3 for time T = 1.2 and = 0.05 in all the modes. The FDM 

offers an improvement in accuracy of several orders of magnitude in the 

error norm over either the MDM or the MAM. The DCM (or in this case the 

FDM of order six) results in over an order-of-magnitude increase in 

accuracy over the FDM (of order four). The advantage of using higher- 

order modal methods (MAM, FDM, or DCM) lies in the ability of those 

methods to approximate the flexibility of the higher, but neglected, modes 

with terms which are functions of the stiffness, mass, and damping 

matrices and the forcing function and, in the case of the FDM and DCM, its 

derivatives with respect to time (see, for example, eqs. (14), (26), and 

(39)). The FDM and DCM offer higher-order approximations by using 

-1 

additional terms together with the pseudo-static response (K Q(T)). For 
this problem FDM assumes N = 4 in eq. (39) and the DCM assumes N = 6 
(hence fourth- and sixth-order, respectively). The moment error norm 
(using five modes) is plotted as a function of time in figure 4. At T = 1.0, 
the value of the moment error norm associated with the MDM is 
equivalent to the MAM value because at time T = 1.0 , Q(T) = 0 and there is 
no difference between MDM (N=0) and MAM (N=l) (comparing eqs. (29) 
and (36)) and, hence, the MDM shows a sharp decrease in error at T = 1.0. 
A similar decrease in error in the MAM qoccurs at time T = 0.8, which 

corresponds to a time when Q (T) = 0. These narrow regions where there is 
a sharp increase in solution accuracy can be anticipated a priori from a 
knowledge of the limes at which the zeroes of the forcing function and its 
derivatives occur. 


A comparison of displacement, moment and shear errors for = 0 and 
times T = 0.6 and T = 1.0 are shown in figures 5a and 5b, respectively. The 
error associated with the displacements is the lowest, the moment errors 
are greater, and the shear errors are the largest. The order of accuracy 
among the displacement, moment, and shear response is expected because 
the moments and shears are functions of successively higher spatial 
derivatives of the displacements (ref. 11). When T = 0.6 and = 0, the 

MAM and FDM (N = 4) are equivalent (fig. 5a) because Q(T) = 0. When 
T = 1.0, the MDM and MAM are equivalent (fig. 5b) for reasons mentioned 
earlier. When T = 1.0 and = 0, the only difference between the MDM and 

• J A a- 2 A‘ 2a 

MAM and the FDM lies in the term ( K MK - <J) Cl il <j) ) Q(T) (eq. (38)). 

The difference between the methods is less at T = 1.0 (Q(T) = 0.) than at T = 

0.6 (fig. 5a), where the difference between the methods lies in the term 
( a a -2aT ^ 

- c|> (J) J Q(T). The higher-order terms are functions of the 

frequencies raised to successively higher negative exponents and, hence, 
should have a negligible effect as higher modes are used providing the 
time-function term does not grow proportionally. 

The effect of damping on the accuracy of the response is shown in 
figures 6a and 6b. Increasing the modal damping does not always 

increase the accuracy of the MAM as suggested in reference 11. For the 

case of a uniformly loaded cantilevered beam subjected to a step loading, 

studied in reference 11, the accuracy of the MAM is enhanced in the 

presence of damping as can be seen from eq. (33). Since all the derivatives 

of the forcing function vanish, the only terms remaining are the pseudo- 

"Cr® rt 

static response and a term which is a function of e . Hence, as C,[ 
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increases, the relative importance of this term on the solution, as compared 

-l 

to the pseudo-static response term (K Q(T)), decreases and, therefore, the 
accuracy of the MAM increases. For other forcing functions, there are 
additional terms which increase in importance as damping increases (see, 
for example, eq. (38) and ref. 13) and so the effect of on accuracy is 

more complex. For example, for the case of a linearly time-varying forcing 
function, the last term in eq. (38) decreases exponentially as the damping 
increases (see eq. (33)); however, the second term increases proportionally 
as increases (the third term vanishes since Q(T) = 0). If the magnitude of 

the second term (eq. (38)) does not decrease with respect to the first term, 
an increase in damping will not necessarily result in an increase in 
accuracy as it does in the unit step function case. Notice that the FDM and 
DCM are much more accurate than either the MDM or MAM for the range 
of time and damping levels considered. 


4.1.12 Step Load.- The step forcing function is discontinuous at time 
T = 0 and hence the integration-by-parts of the convolution integral is 
assumed to begin at time T = 0 + . Including the discontinuity in the 
integration-by-parts results in jump conditions which must be included at 
times during which there are discontinuities in the forcing function and its 
time derivatives (ref. 38). The MAM, FDM, and DCM produce the same 
results for a step forcing function because for T > 0, Q(T) is constant and all 
its derivatives vanish (see, for example eq. (38)). As shown in figure 7, the 
MAM, FDM, or DCM are more accurate than the MDM. A plot of moment 
error, using the first 25 modes, as a function of time is shown in figure 8. 



Over the time range considered (T = 0.001-0.01) the moment errors of the 
higher-order methods are about one-half the magnitude of the error using 
the MDM at time T = 0.001. Also, as time increases, the error associated 
with each method for a given number of modes used tends to decrease and 
at time T = 0.01, the error using the higher-order methods is an order of 
magnitude lower than the MDM. 

For a discontinuous forcing function, such as a unit step, p.(T) at T = 0, 
the MDM exactly predicts a zero response at time T = 0. The higher-order 
methods (N £ 1), however, require a summation of all the modes to predict 
exactly a zero response or the inclusion of appropriate jump conditions. 
Therefore, the MDM will produce qualitatively better results for times near 
T = 0 or close to discontinuities. To calculate the transient response 
accurately at very small times, a large number of modes is necessary. The 
displacement distribution for a unit step loading at time T = 0.0002 and 
C l \ = 0.05 which was calculated using 25 modes is shown in figure 9. As 

expected, the MAM, FDM, and DCM results are equivalent and more 
accurate than results obtained using the MDM. The moment error norm 
associated with the higher-order methods is exceedingly large, at time 
T = 0.0002, when fewer than five modes are used to approximate the 
response (Fig. 10). If a sufficient number of modes is used to predict the 
displacement distribution accurately (m > 24), the higher-order methods 
appear to give better results. 

Hence, for discontinuities in the forcing function and its time 
derivatives, the higher-order modal methods should include the 
appropriate jump conditions. The jump conditions are necessary because 
the integration-by-parts of the convolution integral requires that the 


functions and their derivatives are continuous. If the jump conditions are 
not accounted for in the higher-order modal methods, there will be 
solution errors close to the time of the discontinuity as seen in the 
previous problem. 


4 .1.2 Multispan Beam 


The second problem studied is a simply-supported, uniform multispan 

beam (10 equal-length spans) subject to two loading distributions and one 

/ El 

forcing function. A nominal frequency, o) 0 = '\l ^4 rad/sec, is used to 

normalize time, t, such that T = co () t. The first 30 normalized natural 
frequencies of the beam are listed in Table 2. An analytical solution for 
the mode shapes and frequencies of multispan beams was obtained using 
equations from reference 41. This problem was selected because the 
frequencies are closely spaced (in groups equal to the number of spans, 10 
in this example) and the chances of a neglected higher mode having a 
considerable effect on the response is increased. 


4 .1.2. 1 Uniform Quintic Time-Varying Load.- For this case the load 

distribution is uniform and varys in time as Q(T) = 1000(T 4 -T^), where T is 
the normalized time. The moment distribution, normalized by the 
maximum value of the moment M of the multispan beam at T = 1,2 and 
= 0.05 is shown in figure II. The FDM and DCM are accurate even when 

only one mode is used, whereas the MDM and MAM require 30 and 10 



modes, respectively, for acceptable accuracy (e < 0.01). The spatial 

moment error norm, e, as a function of the number of modes for T = 1.2 is 

shown in figure 12. The first nine modes are nearly orthogonal to the 

T 

uniform load distribution, hence the modal load <J> r Q(T) is negligible and 
has a negligible effect on the MDM response (see fig. 12). The 10th and 
30th modes, however, have an effect on the solution as shown in figure 12 
The effect of these higher modes, however, is taken into account to some 
degree by the pseudo-static response (note the MAM curve for m < 10) and 
to a greater degree by the higher-order approximation of the neglected 
modes used in the FDM and DCM methods. Hence, the FDM and DCM using 
only one mode calculate a more accurate moment response than the MDM 
using 49 modes or the MAM using nine modes. The spatial moment error 
norm is shown in figure 13 as a function of time for each method (using 10 
modes). The accuracy of the FDM is at least two orders of magnitude 
greater than the MDM and at least one order of magnitude greater than the 
MAM. The DCM is, in general, more accurate than the FDM. The MDM and 
MAM are equivalent at T = 1.0, as expected, because Q(T) = 0. At T = 1.0 
the error associated with the MDM decreases an order of magnitude and 
that associated with the MAM increases an order of magnitude as shown. 


4.1 ,2.2 Discrete Quintic Time-Varying Load.- For this loading distribution, 
the solution does not converge in a step-like manner as in figure 12 but 
does so gradually and at a slower rate as shown in figure 14. This 
convergence occurs because the loading distribution is not nearly 
orthogonal to many mode shapes and hence the associated modal load. 
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T 

<j> r Q(T), is not negligible as it was for the uniform distribution. Once again, 
the FDM and DCM converge more rapidly than either of the other, lower- 
order modal methods. The DCM converges faster than the FDM for T > 0.2. 
The DCM and FDM both experience similar error at times T = 0.2 shown in 
figure 15. As time progresses, the DCM appears to be more accurate than 
the FDM. Once again, the accuracy or these higher-order modal methods 
(FDM and DCM) are several orders of magnitude greater than the lower- 
order methods (MDM and MAM). In addition, the accuracy of the higher- 
order methods tends to increase as T increases (see fig. 15). 


4 7 


4.2 Non-Proportional Damping 


4.2.1 Two Degree-of-Freedom Spring-Mass-Damper System.- A simple, 
two-degree-of-freedom spring-mass-damper problem (see fig. 16) with a 
sinusoidal forcing function was analyzed to compare the accuracy of the 
MDM, MAM, FDM, and DCM. For this problem Mj = M 2 = l Kg, spring 
constants K, = K 2 = 1000 N/mm, and damping constant C = 1 N-s/mm. This 
problem was also investigated in reference 39 and included in that 
reference are the particular solutions for polynomial as well as sinusoidal 
forcing functions. The sinusoidal forcing function, sin(Wft), is applied to 

the second mass as shown in figure 16. The natural frequencies are (O^ 
19.54 rad/s and (0 2 = 51.17 rad/s. The system is proportionally damped (A 
is diagonal in eq. (38)) if a = K]/K 2 (ref. 42). Hence, proportional damping 

occurs for a value of Ot = 1.0. The accuracy of each method is assessed by a 
time-integrated error norm (eq. 46) where uj(t) is the calculated response 

using all the modes and uj(t) is the approximate response using a subset of 
the modes. Results were calculated using both the real and damped modes 
(eqs. (26) and (39), respectively). The modal coordinates, q and Z, were 
calculated by numerically integrating equations (7) and (30), respectively, 
using a Runge-Kutta method. For this problem, the FDM used was of order 
four (N = 4 in eqs. (26) and (39)). The time, X, selected for integrating the 
error was chosen to be x = 1 6tc/co f . 
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Results of Ihc error as a function of the forcing frequency using one real 
mode for the proportionally-damped case (a =1.0) is shown in figure 17. 
In general, the forcing function frequency must be lower than the highest 
natural frequency used in the approximate modal response for accurate 
results. As shown in figure 17, the accuracy increases as the order of the 
modal method increases. The results for the lower range of frequencies are 

shown more clearly in figure 18, which is an expanded error scale of figure 

17. The results for the FDM (N=4) and DCM are similar and more accurate 
than the lower order methods such as the MDM (N=0) and the MAM (N=l) 
for C0f< 20 rad/s. For C0f> 20 rad/s, the DCM remains slightly more 

accurate than the FDM; however, as the forcing frequency approaches the 
second natural frequency ( CO f = 51.17 rad/s), all methods produce 

inaccurate results. Results using one real mode or two damped modes are 
identical for the proportionally-damped case. 

Results for the non-proportionally-damped case (a = 20) using two 
damped modes are shown in figure 19. Results are similar to the 

proportionally-damped case with the exception that the DCM exhibits 

surprisingly good results at forcing function frequencies close to the second 
natural frequency. This result is unexplained at present and is believed to 
be fortuitous and, hence, it is recommended that all modal methods should 
include modes whose frequencies exceed the frequency of the forcing 
function. A comparison of the damped-mode solution (using two damped 
modes (eq. (26)) and the undamped solution (using one real mode (eq. 
(39)) is shown in figure 20 for the non-proportionally-damped case (a = 
20). As shown in figure 20, the damped-mode solution using two damped 
modes (dashed lines) produces more accurate results than those using only 


one real mode (solid lines). The damped-mode solution for the FDM and 
DCM are nearly equivalent and result in the smallest error for frequencies 
as large as 30 rad/s. Hence, it may be beneficial, in some cases, to use the 
damped modes to obtain a more accurate solution. 

The FDM produced results that are similar to the DCM results for forcing 
frequencies below the first natural frequency. A comparison of the modal 
methods for a forcing frequency (Of = 10 rad/s is shown in figure 21. Once 

again, the higher-order modal methods result in more accurate solutions. 
The large relative errors near x = 0 are due to the zero initial conditions 
which cause the denominator of eq. (46) to approach zero at X = 0. As 
explained in reference 13, the increase in accuracy with the order of the 
modal method is due to the addition of terms which are functions of the 
generalized stiffness and mass matrices and the force vector and its time 
derivatives. These additional terms approximate the effect of the higher 
modes which were neglected in the modal summation. 


4.2.2 Multispan Beam with Discrete or Uniform Damping and Uniform 
Loading - A multispan beam (five spans) with discrete or uniform damping 
and uniform loading, is shown in figure 22. A similar problem was also 
studied in references 22 and 23. This problem was chosen because it not 
only includes the effect of nonproportional damping, but also the effect of 
discrete damping. The load history was changed in the present study, to a 
uniform quintic varying load. The modulus of elasticity E = 10 x 10 5 
lbs/in. 2 , density p = .28 lbs/in. 3 , and damping constants are Cj = 0.008 lbs- 
sec/in. and Cj = 1.2 Ibs-sec/in. The finite element method was used to 
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discretize the problem in space. Three conventional beam elements (ref. 
23) per span were used to represent the five-span beam as shown in 
figure 22. Results of moment error for the case of discrete damping at 
time t = 1.2 sec are shown in figure 23. Similar to previous results, the 
higher-order methods are consistently more accurate; the FDM and DCM 
have errors e < 0.01 using only one mode, whereas the MAM and MDM 
require 5 and 15 modes, respectively, for comparable accuracies. The DCM 
is shown to be an order of magnitude more accurate than the FDM at time 
T = 1.2. Similar results for a proportionally damped case where = 0.05 
are shown in figure 24. For the case of proportional damping, the DCM and 
FDM results are very similar; for the case of discrete damping (fig. 23), the 
DCM results are almost an order-of-magnitude better than the FDM. 


Table 1, 


Natural 

f 

a>o = 

v 


frequencies of a cantilevered beam 

v rad/sec / 


NORMALIZED 

MODE NATURAL FREQ. 
NUMBER © /(Qq 


1 

3.52 

2 

22.03 

3 

61.70 

4 

120.90 

5 

199.86 

6 

298.56 

7 

416.99 

8 

555.17 

9 

713.17 

10 

890.73 

11 

1088.12 

12 

1305.26 

13 

1542.13 

14 

1798.74 

15 

2075.08 

16 

2371.17 

17 

2687.00 

18 

3022.57 

19 

3377.87 

20 

3752.92 

21 

4147.70 

22 

4562.22 

23 

4996.49 

24 

5450.49 

25 

5924.23 

26 

6417.71 

27 

6930.93 

28 

7463.89 

29 

8016.59 

30 

8589.02 



Tabic 2. 


Natural frequencies of a simply-supported multispan 
beam (10 spans) 


f 

/ El ] 

co 0 = ^ 

V pAL4 rad/sec J 


MODE 

NUMBER 

NORMALIZED 
NATURAL FREQ. 
©/ % 

1 

9.87 

2 

10.15 

3 

10.94 

4 

12.17 

5 

13.70 

6 

15.42 

7 

17.25 

8 

19.07 

9 

20.75 

10 

21.91 

11 

39.48 

12 

40.07 

13 

41.72 

14 

44.11 

15 

46.90 

16 

49.96 

17 

53.12 

18 

56.23 

19 

58.95 

20 

60.95 

21 

88.83 

22 

89.78 

23 

92.23 

24 

95.74 

25 

99.87 

26 

104.25 

27 

108.88 

28 

113.29 

29 

117.13 

30 

119.97 
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Displacement distribution, 


Modes 


Method 
o MDM 

n MAM 

O FDM 

A DCM 


e 

1 0.2890 

1 0.0407 

1 0.0008 

1 0.0011 



a) Displacement distribution 

Figure 1. Representation of displacement and moment errors in 
a cantilevered beam with tip loading using spatial 
error norm e (Q(T) = lOOOfT* - T5) lbs., T = 0.4, and 
damping ratios ^ = 0.05). 
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Normalized moment distribution, M/M root 



Method 

Modes 

e 

o 

MDM 

2 

0.3950 

□ 

MAM 

1 

0.1190 

O 

FDM 

1 

0.0023 

A 

DCM 

1 

0.0033 


- MDM 

30 




b) Moment distribution 
Figure 1. Concluded. 
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Load history 
Q(T), lbs 



a) Force history 

Cantilevered beam with zero damping (£ 
subject to a tip load Q(T) = 1000(T 4 - T 5 ) lbs. 


Figure 2. 



np displacement, in. 


-5 


Q(T) = 1000 (T - T D ) lbs 

I L 1 



b) Tip displat 


Figure 2. 




Moment error norm, 







Moment error norm, 



Figure 4. Variation of moment errors as a function of time, using 
5 modes in the modal summation. Cantilevered beam 
with a tip load Q(T) = 1000(T* - T 5 ) lbs. and where 
damping ratios = Q.05. 



Error norm, 








Shear error norm, 









Moment error norm, 






Moment error norm, 



Figure 8. Variation of moment errors as a function of time, 
assuming 25 modes are used in the modal summation, 
for a cantilevered beam with a unit step tip loading 
(Q(T) = Ll(T) lb. and damping ratios = 0.05). 


Displacement, 


o Ml 

□ M/ 

O FC 

A D( 



x/L 

Figure 9. Displacement distribution in a uniform cantilevered 
beam with a unit step tip loading at time T = 0.0002 
(Q(T) = (1(T) lb. and damping ratios = 0.05). 
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Number of modes 


Figure 10. Comparison of displacement errors for a uniform 
cantilevered beam subject to a unit step tip loading at 
time T = 0.0002 (Q(T) = \i(T) lb. and damping ratios 

= 0.05). 




Q(T) = 1000 (T 4 -T 5 ) lbs 


T = 1.2 


hH 10 spans 

l-*-x 



Figure II. Comparison of normalized moment distributions along 
a simply-supported multispan beam (ten equally- 
spaced spans) and a uniformly-distributed load 
varying with time as Q(T) = 1000(T 4 - T 5 ) lbs./in., 
where T = 1.2 and damping ratios £i = 0.05. 
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Moment error norm 



Figure 12. Comparison of moment errors of a simply-supported 
r multispan beam (ten equally-spaced spans) subject to 

a uniformly-distributed load varying with time as 
Q(T) = 1000(T 4 - T 5 ) lbs./in., where T = 1.2 and 

damping ratios = 0.05. 


Moment error norm, 



Figure 13. Variations of moment errors as a function of time, 
using ten modes in the modal summation, for a 
simply-supported multispan beam (ten equally-spaced 
spans) subject to a uniformly-distributed load varying 
with time as Q(T) = 1000(T 4 - T 5 ) lbs./in. (where 
damping ratios £} = 0.05). 
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Moment error norm, 



Figure 14. Comparison of moment errors of a simply-supported 
multispan beam (ten equally-spaced spans) subject to 
two concentrated loads spaced about the center of the 
first span and varying in time as Q(T) = 1000(T 4 - T ) 
lbs., where T = 0.4 and damping ratios £ j = 0.05. 


Moment error norm 



Figure 15. Variation of moment errors as a function of time, using 
10 inodes in the modal summation, for simply- 
supported multispan beam (ten equally-spaced spans) 
subject to two concentrated loads spaced about the 
center of the first span and varying in time as Q(T) = 
1000(T 4 - T 5 ) lbs., where damping ratios = 0.05. 



M 1= M 2 = 1 Kg 
K 1 = K 2 = 1 000 N/mm 

C = 1 N-sec/mm 
a = 1 - Proportional damping 
a 4 1 - Non-proportional damping 


Figure 16. 


Two-degree-of-freedom spring-mass-damper system. 





Modal Methods (1 real mode 
or 2 damped modes) 



10 20 30 40 50 


a>f , rad/sec 

Comparison of various modal methods as a function of 
forcing frequency for a proportionally-damped two- 
degree-of-freedom problem where a = 1. 


Modal Methods (1 real mode 
or 2 damped modes) 

Mode displacement (MDM) 

Mode acceleration (MAM) 

Force derivative (FDM) 


(4th order) 



to rad/sec 


Figure 18. Comparison of various modal methods as a function o 
forcing frequency for a proportionally-damped two- 
degree-of-freedom problem where a = 1. 



Modal Methods (2 damped modes) 


Mode displacement (MDM) 



tOj , rad/sec 

Figure 19. Comparison of various modal methods as a function of 
forcing frequency for a non-proportionally-damped 
two-degree-of-freedom problem where a = 20. 
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Approx, soln. using 

1 real mode 

— Approx, soln. using 

2 damped modes 



to f , rad/sec 

Figure 20. Comparison of damped- and natural-mode solutions 
for a non-proportionally-damped two-degree-of- 
freedom problem where Ot = 20. 
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t (sec) 

Figure 2 1 . Comparison of various modal methods as a function of 
time for a non-proportionally-damped two-degree-of- 
freedom problem where (X = 20. and (Of = 10. rad/sec. 
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K- 2 . 0 " ■di 

h = .0625" 

T 


Q (t) = 1 000 (t 4 - t 5 )lbs/in. 


Beam cross section 


\ \ 1 1 1 1 1 H H I H I 


3 finite-elements per span 


5-span beam 


E = 10 x 10 5 lbs/in. 2 
p = .28 lbs/in. 3 
C-j = .008 sec - lbs/in. 


C 2 = 1 -2 sec - lbs 
Material properties 


or = 0.05 


Figure 22. 


Five-span beam with applied uniform, quintic time- 
varying load, Q(T) = 1000(T 4 - T 5 ) lbs./in„ and discrete 

dampers. 
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Moment error norm, 



0 2 4 6 8 10 12 


Number of modes 

Figure 23. Comparison of moment errors of a simply-supported 
multispan beam (five equally-spaced spans) subject to 
a uniformly-distributed load varying with time as 
Q(T) = 1000(T 4 - T 5 ) lbs./in., at T = 1.2 and with 
discrete damping. 



Moment error norm, 


Modal method 

— MDM 

- MAM 


FDM 



Number of basis vectors 


Figure 24. Comparison of moment errors of a simply-supported 
multispan beam (five equally-spaced spans) subject to 
a uniformly-distributed load varying with time as 
Q(T) = 1000(T 4 - T 5 ) lbs./in., where T = 1.2 and 

proportional damping with damping ratios Ci = 005 • 


Chapter 5 


Thermal Analysis 


Previous attempts at using mode superposition methods to solve 
transient, linear thermal problems (refs. 31 to 35) have been unsuccessful. 
Thermal problems typically exhibit a wide spectrum response where very 
high frequencies are excited and, hence, a prohibitively large number of 
"thermal modes" are necessary for an accurate solution. Results of 
references 32 and 33 indicate that Lanczos vectors can be effective 
reduced basis vectors for solving linear and non-linear transient thermal 
problems. Since the accuracy of the Lanczos vectors is comparable to that 
of the MAM for structural dynamic problems^ it was expected that higher- 
order methods, such as the MAM, FDM, and DCM, would be effective in 
solving complex thermal problems also. 

Higher-order modal methods, developed in Chapter 3 are used to solve a 
simple linear, transient heat transfer problem of a rod heated at one end. 


The higher-order methods will be shown to be effective in significantly 
reducing the number of modes necessary to represent an accurate 
response. 


5.1 Rod Heated at One End 


The thermal problem selected to study is similar to that presented in 
reference 32. The rod is heated at one end and the temperature at the 
opposite end is constrained to zero (see fig. 25). The forcing function is a 
ramp heat load at one end which ramps up from zero to a peak value at 
time t = 10 sec and down to zero at time t = 20 sec as is shown in figure 
25. The value of temperature at the unheated end is constrained to 0. The 
spacial error norm (eq. 45) is used to evaluate each of the modal methods. 

A total of twenty equally-spaced finite elements were used to model 
the problem. Temperature distributions in the rod, calculated using the 
MDM, are shown in figure 26a. The exact solution, using all 20 degrees-of- 
freedom or modes, is illustrated by the solid line. At time t — 10 sec, the 
peak value of temperature is 400 at x/L = 0. An approximate solution 
using only 5 thermal modes underpredicts the maximum value by 50 
percent (T max = 200) and results in unrealistic oscillations in the 

temperature distribution. If 10 and 15 modes are used, maximum 
temperatures are underpredicted by 25 and 12 percent, respectively. This 
error illustrates the inadequacy of the MDM in accurately predicting 
transient temperatures using a reduced set of thermal modes . 


d-2_ 


Temperature distributions, calculated using the MAM, are shown in figure 
26b. The MAM overpredicts maximum temperatures by over 25 percent 

when four modes are used in the solution. Once again, the temperature 
distribution oscillates, as shown by the dashed line, when an insufficient 

number of modes are used in the solution. The MAM requires 

approximately eight modes for a reasonable solution to this problem. The 
FDM also displays oscillations in the temperature distribution when an 

insufficient number of modes is used in the solution, however, results 

converge to an accurate solution using only five modes (fig. 26c). The 

effectiveness of using higher-order modal methods for reducing the size 
and computational effort of thermal problems is illustrated in figure 27. 
The FDM or DCM require about 28 percent of the number of modes as 
compared to the MDM and about 63 percent of the number of modes as the 
MAM for an accurate thermal response. The ability of the higher-order 

modal methods to predict the transient thermal response accurately using 

very few degrees of freedom or modes, highlights the potential usefulness 
of these methods in reducing the computational size and effort required to 
solve transient thermal problems. 
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Spatial error norm 
e = 

where T a =» approximate temperature vector 

Figure 25. One-dimensional heat conduction problem: rod heated 

at one end. 
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Time = 1 0 sec 



c) Force derivative method (FDM) and dynamic 
correction method (DCM) 

Figure 26. Concluded. 
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T a = approx, temp, vector 



Number of modes 

Figure 27. Convergence of modal methods for a one-dimensional, 
transient thermal problem at t = 10. sec. 




Chapter 6 


Concluding Remarks 


A unified means for creating higher-order modal methods is presented 
which results in algorithms that have increased accuracy when solving 
transient, linear, thermal and structural problems. This new method is 
called the force-derivative method (FDM) because it is based on terms 
which are functions of the forcing function and its time derivatives. 
Several variations are presented for deriving the FDM for both a first- 
order system of equations and a second-order system. One representation 
results in expressions for either a first- or second-order system, which are 
well suited for inclusion into existing finite-element computer programs. 
The FDM unifies previously presented, lower-order methods such as the 
familiar mode-displacement method (MDM) and mode-acceleration method 
(MAM) of structural dynamics and is shown to be equivalent to another 
mode summation method, the dynamic-correction method (DCM), under 
certain conditions. In addition to structural dynamics, it is shown that the 
higher-order modal methods can be used to effectively reduce the size of 
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transient, linear thermal problems, where previous attempts, using the 
MDM, were fruitless. Approximate methods, such as the reduced basis 
methods, which are highly desirable for several reasons: they can 

drastically reduce the problem size and, hence, storage requirements and 
enable the solution of very large problems; and they can decrease the 
computational time and enable large optimization problems to become 

tractable. 

These newly developed higher-order methods have been evaluated 
solely on the basis of their ability to decrease the number of degrees-of- 
freedom necessary to achieve a desired or predetermined accuracy. 
Accuracy has been determined quantitatively, using spatial and time- 
integrated error norms. The accuracy and convergence histories of various 
modal methods are compared for both uncoupled (proportionally-damped) 
structural problems as well as thermal problems and coupled (non- 
proportionally-damped) structural problems. The example problems 
which assume proportional damping include: a cantilevered beam 

subjected to a quintic, time-varying tip load and a step tip load and a 
multispan beam subjected to both uniform and discretely-applied quintic 
time-varying loads. Examples of problems with non-proportional damping 
include: a simple two-degree-of-freedom spring-mass system with discrete 
viscous dampers subjected to a sinusoidally-varying load and a multispan 
beam with discrete viscous dampers subjected to a uniformily-distributed, 
quintic time-varying load. The thermal example problem is a rod 
subjected to a linearly-varying tip heat load at one end with a constrained 
constant temperature at the opposite end. 
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The example problems were chosen to evaluate the effect of: 1) 

spatially and temporally discontinuous forcing functions, 2) forcing 
functions with vanishing higher derivatives, 3) discrete (non-proportional) 
damping and uniform proportional damping in all modes, 4) using 
"damped" or first-order eigenmodes compared to using "undamped" or 
second-order formulation, 5) the level of viscous damping, 6) closely- 
spaced frequencies (which occur in the case of repeated structural 
elements such as the multispan beam example), and 7) using higher-order 
modal methods in solving transient thermal problems. 

In general, the FDM was found to be more accurate than either the MDM 
(zeroth-order method) or the MAM (first-order method) and results in a 

converged solution using fewer modes. The FDM, assuming an order of 

four, was less accurate than the DCM which is shown to be equivalent to 
the FDM under certain circumstance. The DCM, however, assumes a 
particular solution for the transient portion of the solution exists. For 

problems in which there are a large number of closely-spaced frequencies 
(e.g., large truss-type structures and multispan beams), the FDM is very 
effective in representing the effect of the important, but otherwise 
neglected, higher modes. Results for a multispan beam (10 equally-spaced 
spans) and a uniform, quintic time-varying load indicate the FDM and DCM 
methods both produce accurate results using only one mode as opposed to 
the MDM which required 49 modes and the MAM which required nine 
modes for similar values of error. The MDM and MAM results converge in 
a step-like manner with very little increase in accuracy as one to nine 

modes are used in the solution. This step-like convergence occurs because 
the first nine modes are orthogonal or nearly orthogonal to the uniform 

9 2 



loading distribution and hence produce a negligible modal load. For the 
case of discretely applied loads, the convergence is not step-like; however, 
the FDM and DCM methods still converge using many fewer modes than 
the MDM or MAM methods. More modes are required for convergence of 
the shear forces in the beam than for moment convergence and more 
modes are required for accurate moment predictions as compared to 

displacements. This is to be expected since the stresses are functions of 
the spatial derivatives of the displacements and the process of 

differentiation tends to magnify errors already existing in the 
displacement calculations. 

At response times close to discontinuities in the forcing function and/or 
its derivatives, the MDM gives qualitatively better results than the higher- 
order methods, when few modes are used in the approximation. However 
when a sufficient number of modes are used to represent the 

discontinuity, the higher-order methods produce accurate results using 
fewer degrees of freedom. Implementation of the FDM or other higher- 

order methods requires the inclusion of appropriate jump conditions at the 

times when discontinuities occur. It was also found that increasing the 

modal damping levels does not always increase the convergence rate of the 
MAM or other higher-order methods as was previously thought. 

A first-order or "damped mode" solution was found to be effective in 
solving a non-proportionally damped two-degree-of-freedom problem. 
Results of this two-degree-of-freedom spring-mass-damper system subject 
to a sinusoidal forcing function indicate that, for the proportionally 

damped case, a solution using the first-order or "damped" form of the 

equations and two damped modes produces identical results as the second- 



order or "undamped" form using one natural mode. Hence, there is no 

advantage in using a damped or first-order form to solve a proportionally- 
damped problem. However, for the non-proportionally damped problem, 
the use of two damped modes produces more accurate results than using 
only one mode and the second-order form. The DCM has the lowest 
percentage error of all the mode-superposition methods over the 
frequency range of 2 to 50 rad/sec. The FDM produced similar results to 
the DCM up to a forcing frequency of about 35-40 rad/sec. For the 
proportionally damped problem, all the methods were inaccurate near a 
frequency of 50 rad/sec (close to the second natural frequency of the 
system). A multispan beam with discrete dampers (non-proportionally 
damped) was studied and the results of the study indicate that the FDM 
and DCM converge using only one mode as compared to 15 for the MDM 
and 5 for the MAM. 

The higher-order modal methods, such as the FDM, were found to be 
very effective in solving a simple one-dimensional thermal problem of a 
rod heated at one end. Previously, modal methods have been inefficient in 
solving thermal problems because the nature of the problem requires the 
inclusion of almost all the modes for an accurate solution. The ability of 
the FDM to approximate the effects of the higher, but neglected, modes 
results in an accurate solution using only five modes, out of a total of 
twenty modes for the rod problem, as compared to the MDM which 
required 18 modes and the MAM which required eight modes for an 

accurate solution. 
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